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ABSTRACT 


A  set  of  Hermitian  equations  is  constructed  using  a  mode 
description  of  the  dynamic  equations  of  the  ocean  surface. 

These  equations  are  sufficiently  flexible  to  include  the  coup¬ 
ling  of  the  ocean  surface  to  the  wind,  viscous  damping  and 
the  effects  of  surface  tension  as  well  as  the  non-linear  inter¬ 
actions  between  surface  waves.  The  system  of  equations  is  ex¬ 
act  and  from  them  a  system  of  approximate,  first-order  (in  time) 
finite  differential  equations  is  derived  and  solved  numerically 
The  solutions  to  these  coupled  equations  p-  ovide  one  with  a  de¬ 
tailed  view  of  the  growth  of  the  non-linear  surface  waves  and 
energy  spectrum  with  time. 
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1 .  INTRODUCTION 


The  problem  we  wish  to  discuss  in  this  report  is  how  to 
describe  the  physical  structure  of  the  ocean's  surface  by 
means  of  a  direct  calculation  of  the  non-linear  interactions 
between  gravity  waves.  Some  of  the  interactive  mechanisms  them¬ 
selves  had  not  been  identified  until  the  last  decade;  so  it  is 
no  surprise  that  there  does  not  exist  as  yet  a  cohesive  theory 
which  presumes  to  calculate  the  surface  structure  of  the  ocean 
from  first  principles.  The  extreme  difficulty  of  the  problem 
xies  in  the  nonlinear  character  of  the  interaction  process 
which  necessitates  the  construction  and  solution  of  non-linear 
models. 

Theoretical  models  describing  the  ocean  surface  and  the  interac 
tion  of  surface  waves  with  surface  waves,  and  surface  waves  with  in¬ 
ternal  waves,  fall  into  two  rather  broad  categories.  The  first  cate 
gory  comes  under  the  general  heading  of  Ray  Theory,  which  is  space 

oriented  and  primarily  concerned  with  wave  packets,  e.g.,  Whitham 

» 

(1966).  In  this  type  of  theory  one  has  position  and  time  dependent 
wavenumbers  and  frequencies.  The  wave  packet  models  are  also  known 
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as  WKB  or  eikonal  theories.  They  consider  the  ocean  surface 
as  essentially  a  superposition  of  a  number  of  spatially  lo¬ 
calized  wave  packets,  each  distinguished  by  a  characteristic 
wavenumber  k  and  frequency  gj  which  are  related  by  means  of  a 
disperison  relation.  These  wave  packets  move  along  trajectories 
defined  by  the  wave  conservation  equations 


and 


dx 

dt 


V.„  oj 
k 


(1.1) 


(1.2) 


which  defines  the  paths  along  which  energy  naturally  propagates. 

The  second  approach  and  the  one  we  will  use  in  this  paper 
is  a  modal  description  -  the  interaction  process.  Tne  mode 
oriented  models  describe  the  ocean  surface  as  a  superposition 
of  waves  (in  the  inperturbed  ocean  this  would  be  a  superpo¬ 
sition  of  sine  waves) .  Such  models  generally  concentrate 
their  attention  on  the  transfer  of  energy  between  the  different 
modes.  The  most  complete  theory  using  such  a  method  is  that 
due  to  Hasselman  (1961)  .  He  introduces  the  mode  expansion  for 
the  surface  elevation  and  velocity  potential  in  the  dynamic 
equation  for  the  ocean  surface,  as  is  shown  in  Section  2.  In 
acMition,  each  mode  is  expanded  in  a  perturbation  series  using 
the  slope  of  the  waves  as  a  smallness  parameter.  The  non¬ 
linear  interactions  are  given  at  each  order  in  perturbation  theory 
by  the  product  of  an  appropriate  number  of  first  order  terms.  By 


) 


assuming  that  the  sea  state  is  initially  Gaussian,  an  expres¬ 
sion  for  the  redistribution  of  energy  between  th^  interacting 
modes  is  obtained. 

In  a  later  paper  (1966) ,  Hasselman  structures  the  modal 
problem  so  as  to  make  applicable  the  methods  of  Field  Theory. 
The  non-linear  interactions  enter  this  structure  as  pertur¬ 
bation  diagrams  of  increasing  order,  just  as  in  Quantum  Field 
Theory  and  Nuclear  Physics.  However,  all  the  difficulties  of 
convergence  and  mathematical  opacity  also  attend  this  approach, 
making  practical  calculations  without  the  assumption  as  to  the 
statistical  nature  of  the  sea  state  uncertain. 

In  the  limiting  case  of  small  amplitude  waves,  both  treat¬ 
ments,  i.e.,  rays  and  modes,  are  equivalent  since  the  require¬ 
ments  for  validity  of  the  WKB  approximation  are  usually  well 
satisfied  for  the  waves  of  interest.  This  equivalence  is  often 
obscured  in  the  development  of  models  since  spatial  structure 
is  bound  up  in  the  details  of  the  phase  relationship  between 
different  wavenumbers  in  a  mode  description;  whereas  the  wave- 
number  structure  is  dependent  on  the  details  of  the  spatial 
correlations  function  in  the  wave  packet  description. 

In  Section  2  our  model  ocean  is  defined  in  terms  of  Ber¬ 
noulli's  equation  for  the  ocean  surface  and  the  kinematic  boun¬ 
dary  condition  at  that  surface.  Ir.  a  mode  description  these 
equations  are  reduced  to  a  set  of  coupled  equations  which  have 
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only  first  order  derivatives  in  time.  The  method  of  keeping 
the  first  order  simplicity  of  the  equations  is  to  introduce  a 
set  of  hermitian  variables.  These  variables  essentially  de¬ 
couple  the  equations  into  a  set  of  right  and  left  moving  waves, 
with  the  amplitudes  q^+^  and  q^  \  respectively.  The  rapid  time 
oscillating  part  of  the  mode  (Fourier)  amplitudes  for  the  vel¬ 
ocity  potential  (0^)  and  surface  elevation  (C^)  is  in  this  way 
removed  so  that  q^  ^  and  q^  ^  are  slowly  varying  functions  of 
time.  These  new  mode  amplitudes  refer  to  a  Fourier  expansion 
in  both  space  and  time. 

The  utility  of  such  a  representation  of  the  surface  ele¬ 
vation  and  velocity  potential  is  explored  in  Section  3,  where 
the  original  coupled  equations  are  cast  in  the  form  of  an  un¬ 
coupled  eigenmode  equation.  In  this  representation  the  only 
variation  in  the  q^*^  amplitudes  comes  directly  from  the  non¬ 
linear  interactions  which  are  in  terms  of  multiple  products 
of  these  amplitudes.  For  calculational  expediency  the  exact 
equations  are  expanded  and  the  terms  grouped  so  as  to  give 
second,  third  and  higher  order  interactions.  The  expansion 
parameter  is  again  a  wave  slope,  but  we  expand  only  exponen¬ 
tial  terms  and  do  not  use  perturbation  theory.  This  we  feel 
obviates  some  of  the  difficulty  encountered  in  Hasselman's 
treatment  of  the  problem. 
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In  Section  3B  the  specialized  problem  of  the  interaction 
of  only  surface  gravity  waves  is  discussed.  A  resonance  (fre¬ 
quency  matching)  condition  simplifies  the  eigenmode  equations 


so  that  only  the  third  order  interaction  terms  contribute  to 
the  rate  equations  for  the  system.  A  discussion  of  the  inte¬ 
gration  procedure  to  be  used  in  obtaining  the  numerical  solu¬ 
tion  to  the  problem  is  presented  in  the  appendix.  The  integra¬ 
tion  technique  employed  allows  one  to  take  large  time  steps  on 
a  scale  corresponding  to  the  oscillating  part  of  the  solution. 
This  is  done  by  integrating  the  rapidly  oscillating  part  of  the 
interaction  coefficient  analytically  and  numerically  integrating 
only  the  slowly  varying  part.  This  technique  was  used  sucess- 
fully  [Cohen,  et  a_l  (1971)]  in  the  treatment  of  the  laser  heat¬ 
ing  of  plasmas. 

Section  4  explores  some  of  the  specialized  types  of  in¬ 
teraction,  such  as  self-interaction  and  the  scattering  of  waves 
which  do  not  change  in  wavenumber.  The  connection  with  the  non¬ 
linear  Stokes  wave  is  also  made  here.  The  form  of  the  surface 
contours  and  the  representation  of  the  sea  spectrum  in  terms 
of  the  q^-)-  amplitudes  is  discussed  in  Section  5. 

A  preliminary  calculation  is  presented  in  Section  6,  which 
has  direct  bearing  on  the  experiment  of  Benjamin  and  Feir  (1967) 
in  which  they  determined  a  single  gravity  wave  to  be  unstable. 

A  single  mode  of  large  amplitude  is  allowed  to  interact  with 
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its  side  bands  which  are  an  order  of  magnitude  smaller  in  am¬ 
plitude.  It  is  observed  in  the  calculation  that  energy  dif¬ 
fuses  out  of  the  primary  mode  into  the  other  modes  in  the  sys¬ 
tem,  until  the  neighboring  and  primary  modes  are  of  comparable 
amplitude.  Energy  is  then  transferred  rapidly  between  all  the 
modes  in  this  region  of  k-space.  The  spatial  picture  of  the 
surface  when  this  rapid  transfer  process  takes  place  is  a  break¬ 
up  of  the  primary  wave  into  a  number  of  wave  packets.  This 
breakup  into  packets  is  the  instability  found  by  Benjamin  and 
Feir.  A  detailed  comparison  with  these  results  is  made  in 
the  text. 
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2.  MODE  COUPLING  ANALYSIS 


In  the  following  analysis,  the  ocean  is  assumed  to 
be  homogeneous  and  irrotational ,  that  is,  if  $  is  the  potential 
describing  th.j  velocity  field  (u  =  V$),  then  V*u  =  0.  The 
velocity  potential  also  satisfies  Laplace's  equation 

V2  4>  =  0  ,  (2.1) 


since  the  fluid  is  assumed  to  be  incompressible.  If  we  define 
the  quantities:  p  as  the  fluid  density,  g  as  the  gravitational 
acceleration,  and  p  as  the  pressure,  we  may  write  the  momentum 
equation  of  the  fluid  as 


du _ 1 

3tf  p 


A  P 


+  g 


(2.2) 


In  terms  of  the  velocity  potential,  Equations  (2.2)  may  be 
written  as. 


du  _ 
3t  " 


y$.y$ 


(2.3) 


where  we  have  made  use  of  the  fact  that  the  d/dt  is  the 
Eulerian  derivative 


d 

at 


t  V4>*7 


(2.4) 


the  density  p  is  constant  (assumed)  and  z  is  the  vertical 
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coordinate. 


We  may  integrate  Equation  (2.3)  immediately  to  obtain 

$t  +  jV<J>*V4>+^E+g:p*0  at  z  =  5  (2.5) 

with  the  condition  that  hydrostatic  equilibrium  must  prevail 
at  infinity,  where  Ap  is  the  incremental  pressure  with  respect 
to  ambient  and  z  *=  z,  is  the  free  surface  of  the  ocean.  A 

second  equation  may  be  obtained  by  recalling  that  the  rate  of 
increase  in  the  wave  height  following  a  fluid  element  is  the 
vertical  component  of  the  fluid  velocity: 

It  *  ft  at  z  =  5  •  (2'6) 

Using  the  notation  (V)  for  the  horizontal  gradient,  we  may 
write  Eq.  (2.5)  and  (2.6)  as 

+  7  (W2  +  9C  +Ap/p=  0  at  z  =  t  (2.7a) 

and 

C  +  Vt’V;  -  *  -  0  at  z  »  (  ,  (2.7b) 

t  2 

as  the  set  of  coupled  equations  describing  the  ocean  surface. 

To  obtain  an  expression  for  the  pressure  in  Eq.  (2.7a) - 

we  must  consider  the  fact  that  the  interface  between  two  fluids 

is  in  a  state  of  uniform  tension.  The  pressure  condition  at 

this  interface  can  be  obtained  by  considering  the  vertical 

forces  acting  on  a  strip  of  surface  of  width  <5x.  If  p  is  the 

pressure  just  below  the  water  surface,  pa  the  pressure  of  the 

£1 
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air  and  the  surface  tension,  then 


(P  - 


Pa)«X 


+  6 


T  — 
L1  3x 


(2.8) 


is  the  torce  balance  acting  on  the  strip  6x.  We  therefore 
have 


The  quantity  pa  is  variable  if  there  is  a  wind  blowing  over 
the  ocean  surface.  We  model  this  variation  in  pressure  by 
the  expression 


p"  ^a  "  u  ?x  ^ ambient  '  (2.10) 

which  gives  the  in-phase  pressure  variation  at  the  surface 
due  to  the  action  of  the  wind?  Substituting  Eqs.  (2.9)  and 
(2.10)  into  (2.7a)  yields 


*t  + 


1 

7 


(V4>)2  +  gC  + 


“H 


-  Y 


ill 

3x2 


(2.11) 


where  we  have  assumed  conditions  to  be  completely  uniform  in 
the  y-direction.  i.e.,  a  series  of  ridges. 

In  the  above  equations  (2.11)  and  C2.7b),  we  have  assumed 
the  existence  of  a  velocity  potential  (u  *  V$)  based  on  the 
assumption  that  the  ocean  is  irrotational .  We  have  also  pro- 


*  Equation  (2.10)  is  intended  to  indicate  how  one  includes  the 
contribution  from  the  wind  in  the  dynamic  equations,  the  de¬ 
tailed  form  of  the  term  will  change  in  a  realistic  calcula¬ 
tion,  e.g.,  making  a  a  function  of  k  as  in  the  Miles-Phillips 
Model. 
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vided  a  means  by  which  energy  is  supplied  to  the  surface  waves, 
i.e.,  by  the  action  of  the  wind.  We  must  now  discuss  the  mech¬ 
anisms  by  which  this  energy  supplied  by  the  wind  may  be  dissi¬ 
pated.  The  most  dramatic  process  is  the  formation  of  "white 
caps"  or  wave  breaking.  In  our  model  this  breaking  directly 
transfers  energy  into  heat  and  therefore  out  of  the  wave  sys¬ 
tem  and  is  not  directly  included  in  the  above  equations.  The 
"white  capping"  should  appear  as  instabilities  in  the  above 
solutions  when  the  wave  slope  has  exceeded  some  critical  value. 

A  second  mechanism  for  the  dissipation  of  energy  is  the  action 
of  molecular  viscosity,  which  has  not  as  yet  been  included  in 
the  above  equations. 

To  model  the  effect  of  viscosity,  we  first  transform  the 
system  into  a  modal  description  by  introducing  the  plane  waves 

Xk(x)  «  elkx//L  (2.12) 

where  L  (*  length)  is  a  dimension  of  interest  tor  the  problem, 
and  decompose  the  vector  potential  and  wave  amplitude  into  their 
Fourier  components, 


and 


*(x,z,t)  »]Cxk(x)  elkl*  4>k(t) 

Jc 


C (x,t) 


<k(t) 


(2.13) 


(2.14) 
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Introducing  Eqs.  (2.13)  and  (2.14)  into  Eqs .  (2.11)  and  (2.7b) 
yields  the  set  of  coupled  equations  in  terms  of  the  modal 
amplitudes, 


pCk(x>  j*k  e!*l2  +  gck  +  iokCk  -  Yk2!R 

=  ( V4> ) 2  ;  z  ~  £ 


(2.15) 


and 


V0*VC;  z  =  ;  .  (2.16) 


We  wish  to  rewrite  Eqs.  (2.15)  and  (2.16)  in  a  form 
where  the  left  hand  sides  of  the  equations  are  linear  and  the 
right  hand  sides  contain  all  the  non-linear  coupxing  effects. 
To  do  this  we  evaluate  (VO)  and  V0*V£  to  be  as  follows: 


(VO)2  =  02  +  02 

x  z 


and 


VO*  Vc 


'l)  v. 


(Ul  +  |p|) 


x*x. 


(2.17) 


Z\  l\z 

*P*£«pe  X4Xp  (2.18) 

l,p 


and  remove  the  exponentials  to  the  right  hand  side  of  Equations 
(2.15)  and  (2.16) .  Then  by  multiplying  these  equations  by 
x£(x)  and  integrating  over  x,  we  obtain 
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+  95 k  =  |  2  <*p  -  1*1  IpI>vp 

£,p  P 

/e^l  +  |p|)5(x)^(x)Xp(x)x*  (x )  dx 

"2  -1)  x£(x)x^(x)dx 


-  iakck  -  vk  ck 


(2.19) 


and 


!■  1 1. 1  a 

”k  '^k 


■  £  *p*t«p  jlu|Ux) 

£  ,p 

+  SlAUA/(e|A  l?(x)-l)x£Xkdx 


X„  X~XiQx 

x,  P  K- 


(2.20) 


or  defining  F1 (k)  and  F2 (k)  by  the  right  hand  sides  of 
Eqs .  (2.19)  and  (2.20),  respectively,  we  have  the  set  of 
equations 


Wk =  Fi<l5) 


where  Gk  =  g  +  yk 


and 

tk  -  |k|4>k  "  F2 (k) 


(2.21) 


(2.22) 


which  to  this  point  provides  a  formally  exact  description 
of  the  interaction  process  in  f  single  horizontal  dimension. 
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We  now  return  to  the  discussion  of  the  viscous  damping 
of  high  frequency  waves.  In  a  viscous  fluid  the  interaction 
between  water  molecules  produces  shearing  effects  when  the 
water  is  in  motion  and  destroys  the  irrotational  assumption 
imposed  above.  In  such  a  system  we  assume  that  we  can  super¬ 
impose  the  rotational  character  of  the  fluid  on  the  irrotational 
so  that 


u  =  V$  +  u' 


(2.23) 


where 


Vxu'  ^  0  . 


The  energy  of  the  fluid  motion  is  gradually  dissipated  by 
these  shearing  effects  and  may  be  modeled  by  determining  the 
rate  at  which  the  energy  is  being  dissipated. 

Following  Phillips'  discussion  of  viscosity,  we  write 
the  rate  of  strain  tensor  as 


3Xj  3x^ 


(2.24) 


where  (i,j  -  1,2)  label  the  components  of  the  coordinate 
system  and  velocity.  Assuming  the  ocean  to  be  an  isotropic, 
incompressible  Newtonian  fluid,  and  y  to  be  the  coefficient 
of  viscosity,  we  may  write  the  frictional  force  per  unit 
volume  in  terms  of  the  derivative  of  the  rate  of  strain  tensor, 
i  .e . 
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since  V*u  =  0.  We  may  also  write  the  rate  of  working  against 
viscous  forcer  as 


3eii 

uifi  =  2“  ui  ST1 


**  2  y 


e 


e 


where  e  =  lyte^)2  =  2yei;. 


(2.26) 


(2.27) 


The  first  term  on  the  right  of  Eq.  (2.26)  represents  the 
viscous  energy  flux;  the  second  term,  the  rate  of  energy 
dissipation  per  unit  volume  by  molecular  viscosity . 

The  rate  of  energy  dissipation  per  unit  area  at  the 
surface  of  a  very  deep  ocean  can  be  written,  using  Eq.  (2.26) 


(2.28) 


The  contribution  from  the  surface,  layer  is  assumed  to  be 
small,  the  primary  effect  coming  from  the  ocean  interior. 

Using  only  the  irrotational  part  of  the  velocity  from  Eq.  (2.23), 
we  may  write  Eq.  (2.28)  up  to  second  order  as 

.o 


•  f  0  fa2* 

E  -J  2y  3x.3x. 
—00  L 


dz 


(2.29) 


We  evaluate  Equation  (2.29)  using  a  single  mode  solution  to 
the  linear  equation,  i.e.. 
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(2.30) 


$  =  c<^k  ekZ  sin  jk(x-ct)J 
so  that  Equation  (2.29)  is  reduced  to 
E  =  -  2y  <j>2k3  c2 

Using  the  fact  that  in  deep  water 

p  _  1  n  .2  2  , 

E  "  2  p  C  k 


(2.31) 


the  attenuation  coefficient  Yv  for  this  wave  can  be  written  as 

Yv  =  2E  =  2vk2  (2.32) 

where  v  =  y/p  is  the  kinematic  viscosity  of  water.  The  energy 
density  o?;  the  wave  field  decreases  as  exp  j-  2Y^t|  and  the 
amplitude  by  exp  Yvt J.  We  may,  therefore,  model  the  effect 
of  viscosity  in  Equation  (2.21)  by  including  the  term 

“2V  k2$k  (2.33) 

in  F1 (k)  . 


The  equations  ior  the  ocean  surface  now  take  the  form 

^k+Gk^k  =  Fl(k)  "  2vk^k  =  Px(k)  (2.34) 

and 

^k  -  lkl^k  =  F2(k)  =  F2(k>  ,  (2.35) 

which  includes  models  for:  (i)  the  generation  of  waves  by 
wind,  (ii)  the  non-linear  interaction  of  surface  waves, 

(iii)  the  effect  of  damping  by  the  dissipation  of  energy 
through  molecular  viscosity,  and  (iv)  the  generation  of  capil¬ 
lary  waves. 
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3.  EIGENMODE  DESCRIPTION 

The  structure  of  Eqs .  (2.34)  and  (2.35)  suggests  that 
it  may  be  possible  to  transform  these  coupled  equations 
to  a  decoupled  representation.  We  can  in  fact  define  the  new 
variables 


and 


4>k  -  */G^1(k) 


Ck  -  i  /|  k  |  i|»2(k) 


(3.1a) 


(3.1b) 


which,  when  substituted  into  Equations  (2.34)  and  (2.35),  yield 

i^k)  -  /GjJ  kj  \p2(k)  =  F^k)  (3.2a) 

^k 

and 


ij2(k)  -  /GjJk|  ^(k)  =  -~=r-  F2(k)  . 


(3.2b) 


If  we  uake  the  sum  and  difference  of  Eqs.  (3.2a)  and  (3.2b), 
we  obtain 

i(*1+i2)  <w  -  Fl(k)  +  -^=  F2(k) 


and 


i<*r*2)  ^  Vk)  -  Vk> 


which  in  terms  of  the  quantities 
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Preceding  page  blank 


(±) 


±  [♦,<«  *  ♦,<»]  k*  ^  J(^) 

k  ■ 


and 


F(-J  (k) 


—  F.(k)  - 

1 


F2<k)u~r 


iki 

!k|3'2 


1 

H 


1 


(3.3) 


(3.4) 


becomes 


ib<+)  -  CiM"  »i+)  - 


F<+) (k) 


(3.5a) 


and 


k  1  +>^JkT  b<-> 


F(_)  (k) 


(3.5b) 


In  terms  of  the  orthonormal  matrices 


C(+)  =  -A  I1 

/2 


and  C(_)  = 


i(-:) 


Eqs.  (3.3)  to  (3.5)  can  be  written 


b.  -  b<+)C(+)  +  b<",C(">  . 

~k  k  ~  k  — 


(3.6) 


(  +  )  (  +  )  (_) 
F  (k)  *  F  (k)  C  '  +  F  (k)  C '  ; 

»v 


(3.7) 


and 


ib.  -  Ab,  =  F  (k) 
~k  ~*>-k  ~ 


(3.8) 


where 


AC<±>  -  ±  /GjkT  £(±) 


(3.9) 
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Using  the  preceeding  equations,  we  may  express  the 
velocity  potential  and  wave  height  in  terms  of  the  eigenmode 
amplitudes  of  the  problem  [b^^l  as  follows 


Both  the  velocity  potential  and  the  wave  height  are  real 
quantities  so  we  may  equate  them  with  their  respective  complex 
conjugates 


In  a  similar  manner 


*k 

* 

*k 
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1 

7 


and  therefore 


)  =  -  i(b 


<+)* 

-k 


(3.12b) 


Adding  and  subtracting  Equation  (3.12)  yields  the 
following  conditions  on  the  amplitudes 


(3.13) 


We  can  remove  the  rapidly  oscillating  part  of  the  expres¬ 
sion  in  our  eigenmode  equations  [Equation  (3.5)]  by  introducing 
the  variables 


e 


T  iu>,  t 
k 


(3.14) 


where  is  assumed  to  be  a  slowly  varying  function  of 

time.  If  we  equate  to  /GjJ  k  |  in  Eq.  (3.5),  then  we  have 
the  new  equations 


F(±)  (k)/T 


(3.15) 


for  our  system. 

(+) 

To  interpret  the  significance  of  the  ,  let  us 
consider  a  single  mode  propagating  to  the  right.  Then, 
according  to  Eq.  (2.14)  the  wave  height  is 

C(x,t)  =  Xk(x)  Ck(t)  +  X_k (x)  C_k(t) 
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this  becomes 


<?k 


and  b^,  ^ 

-k 

i 

- [»r; 

J  as  the 

(  1  \ 

r (+)/ 

(kx-w^t) 

V  2  |k|  ) 

[*k  e 

t  s  ur1 

and 

■  i^+,i 

A 

t 

this  becomes 

5  (f/fc)  '■  -kj[+)|  sin  jkx  -  o)kt  +  eRJ  .  (3.1 

The  dimensionless  quantity  |q^+)|  thus  represents  the 
ratio  of  the  maximum  wave  amplitude  to  X.  We  expect 
0  <  <  1 '  since  at  the  upper  limit  very  strong  coupling 

between  modes  will  occur.  This  notion  is  consistent  with  the 
analysis  done  by  Stokes  on  finite  amplitude  effects  on  gravity 
waves  in  deep  water.  it  was  determined  that  a  maximum  crest 
angle  (120°)  for  gravity  waves  existed,  after  which  the  wave 
becomes  unstable.  This  restriction  is  easily  maintained  in 
terms  of  the  slope  of  the  gravity  wave,  i.e.,  the  ratio  given 
by  Equation  (3.16).  The  quantity  q^+)  is,  therefore,  a  more 
natural  variable  for  the  problem  than  the  wave  height. 
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3A.  The  Approximate  Equations 

Although  Eq.  (3.15)  is  formally  exact,  in  its  present 
form  it  is  not  of  much  use.  It  may  be  put  in  a  more  manageable 
form  for  numerical  calculation  by  making  a  number  of  reasonable 

A 

approximations.  Consider  the  function  F^(k)  , 

.  /*(l*Mpl)c(x> 

Fl(k)  =  2  S  **p  "  lallpl^j^p  Je  xixpxk  dx 

.  f  /  | * I C (x)  \  *  2 

-L**y(e  -  ljX^dx-  iak  Ck-  2vk  $k 

(3.17) 


in  this  expression  the  function  t  multiplied  by  a  wavenumber 

in  the  exponential  is  a  small  quantity  since  the  slope  cannot 

exceed  the  value  1/7.  This  value  is  a  result  of  the  Stor.es 

analysis  mentioned  earlier.  The  wavenumber  in  the  exponential 

is  generally  restricted  to  a  region  of  wavenumber  space  in  which 

the  surface  elevation  is  being  expanded,  i.e.,  1*1  ^  lxmax 

and  I  k  .  -k  I  <<  111,  so  that  I  £  I C  (x)  has  a  critical  value 
1  min  max  1 

which  is  on  the  order  of  2tt/7  Ml.  It  is  therefore  not  unreason¬ 
able  to  expand  the  exponentials  under  the  integral  sign  in  Equa¬ 
tion  (3.17)  and  keep  terms  to  only  third  order  in  the  amplitudes 
when  one  is  not  considering  capillary  waves  or  those  gravity 
waves  close  to  breaking. 

We  perform  this  expansion  in  Equation  (3.17)  and  replace  Z. 
by  its  Fourier  expansion  [Equation  (2.14))  to  obtain, 
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Fl(k)  »  j  53  Up  -  Ii||p|)$  a  -  53  Uli*CD  -  iakc. 

l+p«k  1  p  &+p°k  p 

SZ  Jl 

+£  7  Up  "  Ullp|)l|t|+  IplH^pC,,-  2vlt2*k 

l+p+n«k 

/ L 

-E  J$J-  i  n  c  .  (3.18) 

l+p+n-k  &  P  n 

/L 

The  time  derivatives  of  the  amplitudes  may  be  replaced  in  this 
expression  by  their  value  from  Eq.  (2.19)  to  obtain 

A  , 

F^k)  «-  iak;k  "  2vk  ?k 

+  ^'p.k|?<lp  ‘  |t|lpl)Vp  +  (=!  +  i“l)l*l«t?p 

~PT~  . 

+  2v^tCpj 

+^+n-k  I*  UP  *|l|lp|)  1,11  +  lpl  *  l‘+Pll*t»p«n 
JT.  . 

*  1^-  l(G(  +  ialUt  *  ^'’^♦nJCpCp  .  (3.19) 

In  a  similar  way  we  obtain  the  following  expression  for  (M 

F2  (k)  *  ^  l(i+p  )$,C  +  X)  |  £p  -Lid—  ||jt‘4>  c  t 

>»P«k  1  p  t+p+n»k  I  1  9lVn  * 

/E 

(3.20) 

Our  eigenmode  equations  are,  however,  in  terms  of  q/ and 

|^\  A  ^ 

F  -  (k)  not  $y,,  and  F(k).  We  must  therefore  determine  F'-* 
in  terms  of  qj,^  . 
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We  recall  that 


F(-)(k)  -  j  -- —  F,(k)  i  - 

*  ^gt  1  /J|TF 


b  ;>“l 


which  may  be  written  using  Eqs.  (3.19)  and  (3.20)  and  Eqs .  (3.10) 


and  (3.11)  ; 


=Tr«««k,.  E  | -ip K’-r’K’-n 

A  • 


t+p*-k 


P  P 


•vSUK’-HK'W' 


k  |\I 


J  a. 
j  tpn  | 


♦  E  |“tp»[bt,+,tbi‘,][bpt,+bp',][b"+,‘b"',l 

t+p+n«k 

•  v[>!"  •■.iw' 

*  C..W’ "W'-ril  (1*1) 

.  u.  tSk-1 -“i'1] 


(3.20) 


-  |a2SL  (lp.,l||p|)|  , 

-  -i  Vui  |p|  (g£  +  iot) 


Y(±) 

Ytp;k 


kGt  */|P| 


-  Jojpl  *  1  *-T-  ViTi 


l(l+p)  » 


-  •  -^-p-  (^p-  i  m  ipi)(iii  + 1 p i  -  i*+pi)  Vi n i /2  * 

2 

-J— 1—  l|  |  p  |  |n  |  /  2  (g£  +  ial^  , 


r 

£pn ;  k 


I  P I  I  n 


(3.22) 


The  complete  expression  for  Equation  (3.16)  is  as  follows: 


vk 


i+p-k 


lip  +  Btp  +  Ytp;k ] 


♦  q(‘V*)<.i(V“p)tTa  +  ft  -v<!>  1 

qi  qp  lalp  +  6lp  Ytp;kJ 


+  ^r^p'’0 


-i (w.-o J  ) t 


1  V "  £ 


„  .  o  v(t)  1 

a£p  Ptp_T£p;kJ 


(-)  (  +  )  iK-°,D)t  r 
+  )qp  e  hP-Bip+Y 


-fi  +Y(±)  (  k 

£p  Dlp  Tip;kJ  \  pT 


jzr  La^pn  Btpn  Y£. 


,+w  it  r  ...  n 

[a£pn  +  Btpn  +  Ytpn?k] 


£+P+n«k 


*  q<+>q(+>q(-,e'i(V“p‘“n)t  T  (1)  1 

Hp  4n  L  *pn  Pipn  Y£.pn;kJ 


/ ,%  /  \  # ,  \  -i (ue-w  +w  )t  r 

+  "  Kn'6 


—  K  —  Y 

Ipn  Ipn  1 Ipn 


*>  1 

pn  ;kj 
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+ 


+ 


+ 


+ 


+ 


± 


<r’We 


i  (ui.  -w  -w  )  t 
*<  p  n 


'  .  v 
_  Apn  pApn  TApn 


J 


i(w  -to  -Ho  )t 

(-)(+)(-)  1  p  n  r.  +B  _Y(±) 

qA  qp  qn  LK*£pn+  PApn  Yj,pn 


qA  qP  qn  e 


A  p  n# 


[- 


+  8  -  Y 

Apn  PApn  1  Apn 


J 

*>  ] 

pn;  k  J 


q  A  qp  qn 


i  (w  A +0)  +w  )  t 


A  p  n' 


-a  -  8  +  Y 

Apn  Apn  Apn 


*>  II 

pn;k J| 


i«K  # 1 2ivic  V¥{ 


/? 


(3.23) 


3B.  The  Resonance  Condition 

We  expect  that  in  the  development  of  the  amplitude  of 
the  kfc^  mode  non-resonant  terms,  being  oscillatory,  will  make 
a  zero,  i.e.,  negligible,  net  contribution  for  the  interaction 
of  surface  gravity  waves  alone.  Only  those  combination  of 
frequency  close  to  i.e.,  those  in  resonance,  will  contrib¬ 
ute  substantially  to  qJ[±)  .  We  therefore  rewrite  Equation 
(3.23)  for  a  right  travelling  wave  as 
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;<+>  _  f oik j kj  _  2  \  „<+) 


2u)i 


-  vk  q 


] 


£.+p«k+n 


a£.p-n  +  ^Jlp-n  +  ^£p-n;k 


„(+)„(+)„(-) 

q£  qp  q-n 


-  h- 


-8.  -Y.<+)  J  q.<+)q(_)q(+) 

np  x-np  T£.-np;kJ  4£.  ^-n  4p 


•  [, 
*  JjZS 


-ft  +  <  +  >  1  a(")a(  +  )a(+) 

npA  B-npi,  Y-npi;kl  q-n  qp  qn 


-.1 


ipn 


exp 


ji(ta,k+a)r1‘aJH-up)t[ 


which  we  can  rewrite  in  the  form 


(+)  .  [aklkl  _  ,,1,2!  „(  +  ) 


2  co, 


-  vk 


t 


Y  ir(  „  qf  +  )q^q<-) 
/  ^  dpn; k  4p  M-n 


fc+p*k+n 


where 


£pn;k 


2  (a£p-n  +  ®£p-n  — ^£p-n;k  af-np  +  ^Jl-np  +  ^f-np;k 


-a  »+S  .  -  v  ^  +  ^  . 

-npi,  -np£  -npf; 


'  k  /  \  2Gk  £pn  J 


and 


(3.24) 


(3.25) 
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(3.26) 


ft-  .  =  OJ.+U)  -10--U) 

Jlpnk  k  n  Z  p 


If  we  define  the  functions 


%n(t> 


„(+>„(+>  r„<+n* 
qt  qP  Lq"  J 


u  *  “ipnk  ' 

then  Equation  (3.24)  can  be  written  in  the  form 


(3.27) 


(3.28) 


•(+)  - 
q,  = 

"ak  k  ,  2 

- vk 

L  2a,k 

where 

r  = 

£pn ;  k 

R£pn;k  +  l3:A 

P  = 

wk  |k| 

x.pn;  k 

16Gk  | Hpn | 5 

lF £,pn ; k  Hfcpn(t)  e 


£+p«*k+n 


(3.29) 


^GiGp  (l  £pI  -  *p)  (l  *1  +  Ip  I  -  I  *+p|)  VM" 

^GzGn  (l*n|  +  *p)(l*l  +  I n  I  -  VTpT 

|/cnGp  (|np|  +npj(|n|  +  |p|  -  |n-p|)  ^\T\ 

(|lpn|)2  (g„  n2  +  2Gt  t2)  +  (n2-2"p)M(|4|)?V5A 

V^k  2J.  (ft  -  n  +  p)  1 1 1  (|f|)?| 
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W  -  ii57  j^7i  [»I*P»1: 2  2‘3'"3 

+  2v^^  n4  |pH|2  -2^1*  |pn|2^  . 

m 


(3.30) 


The  numerical  integration  of  Equation  (3.29)  is  discussed 


in  Appendix  A. 


Preceding  page  blank 
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4.  DISCUSSION  OF  EQUATION  (3.29)* 

The  wavenumbers  in  the  interval  0  <  x  <  L  are  of  the  form 

k  =  4r  K  '  P  =  ^7  P  ,  etc.  ,  (4.1) 

where  K,  P,  etc.,  are  integers.  The  wavenumber  matching  con¬ 
dition  in  Equation  (3.24)  using  Equation  (4.1)  is  then  an  in¬ 
teger  matching  condition. 

The  quantity  T  [Equation  (3.30)]  has  the  dimensions  of 
a  frequency  and  is  of  the  same  order  of  magnitude  as  the  mode 
frequency,  that  is, 

T  =  0(wk)  ,  (4.2) 

where  is  a  characteristic  wave  angular  frequency  in  the  wave 
number  region  of  interest.  We  recall  that  the  q^  are  dimen¬ 
sionless  ratios  of  [from  Equation  (3.16)], 

wave  height  2  tt 
wave  length 

When  the  q^'s  are  nearly  unity,  the  non-linear  interaction  has 
a  time  constant  comparable  to  the  wave  period. 

When  the  q^'s  are  small  compared  with  unity,  the  character 
istic  time  constant  o'"  the  non-linear  terms  is  -  |wk  |  a  |  2  j -1 . 

In  this  Section  we  set  G^^g,  i.e.,  neglect  surface  tension 
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When 


l%n]J  >>  [uklal2] 


(4.3) 


the  oscillating  exponential  will  tend  to  "wash  out"  the  non¬ 
linear  coupling. 

If  we  set  [assume  £,  p,  k,  n  £  0) 


k  =  £  -  n 
n  =  p  +  n, 


(4.4) 


we  satisfy  the  wavenumber  matching  condition  k  +  n  =  £  +  p.  The 
frequency  mismatch  (3.26)  is  then 


=  “  /g  +  /p  -  /£-r)  -  /p  +  r)J  , 


(4.5) 


-  p  <  n  <  £■  . 


Next,  we  assume  £  £  p  and  write 


(4.6) 


q  = 


to  obtain 


=  E 


£  ' 


y  = 


»  -q  <  y  <  i» 


-  —  =  f(q,y)  »  1  +  /q  -  /I=y 


-  /q+y 


(4.7) 


(4.8) 


At  the  y-end  points  f  is  positive  and  has  the  value 


f(q,  -q)  =  f (q ,  1) 


2/q 


1+  /q"  +  /1+q 


(4.9) 
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f  has  two,  and  only  two,  zeros.  These  occur  at 


n  =  o  ,  n  =  A  -  p 


(4.10) 


Between  these  zeros  f  is  negative  and  has  the  minimum  value 

(1  -  ^q)2 


-(q)  -  -  -  jrfr 


/q  +  /2  (1  -  q) 


The  quantity  fmin(<3)  is  plotted  in  Figure  1. 


For  the  self-interaction  of  a  single  mode,  we  have 

k  1  kkk; k  2  1  ^  8g  4o>k  ) J  ‘ 


(4.12) 


For  scattering  with  no  change  of  wavenumber,  we  have 


r  =  r  +  r 

nk  nkn;k  Knn;k 


■  '  xj1  +VM(2  +  n  - 

lt0k  / ok  _  vk2  \(k\2 
2  y  4g  "  2wk  J  \n  / 


1  -  * 
n 


(4.13) 


The  appropriate  form  of  Equation  (3.24)  in  the  case  of 
pure  resonant  interactions  is 


•  (+)  -  l 

%  -  V 


ak  k  ,2 

■= -  -  vk 

2u). 


)«k+ 


+  ir. 


_  ( t ) 
qk 


„(  +  ) 

qk 


♦  Li 

n  (^k) 


i  r  „  q(+> 

nk 


2  <+) 


(4.14) 
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The  quantity  [-  r^/w^]  is  plotted  in  Fig.  (2)  as  a  function 
of  the  ratio  (k/n)  for  the  case  a  =  v  =  0. 


Let  us  consider  the  rate  equation  for  the  self-interaction 
of  a  single  mode  in  the  absence  of  wind  and  neglecting  vis¬ 
cosity,  i.e.  , 


q‘+!  -  4+) 


(4.15) 


where  we  have  used  Equations  (4.12)  and  (4.14).  We  may  pre- 


multiply  Equation  (4.15)  by  q^+^  and  premultiply  the  complex 

.(+> 


conjugate  of  Equation  (4.15)  by  q^  and  add  the  two  equations 
to  obtain 


at  l^+!  I  -  ° 


(4.16) 


which  implies  that  the  modulus  of  q^+)  is  constant.  We  there¬ 
fore  write 


(+)  iXk 

qO  c,  e 


(4.17) 


where  ck  is  constant.  Substituting  Equation  (4.17)  into  Equa¬ 
tion  (4.15)  we  obtain 


w. 


~2  ,ck 


(4.18) 


which  yields 


.(  +  )  _ 


(4.19) 


as  the  solution  to  the  interaction  equation  for  a  single  mode 


[Equation  (4.15)]. 


*  Note  that  this  result  will  not  change  if  the  T  ,  terms  in 

nk 


Equation  (4.14)  are  also  included  in  Equation  (4.15). 
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In  this  single  mode  case  we  can  write  the  velocity  poten¬ 
tial  in  terms  of  its  Fourier  amplitudes  [Equation  (4.19)]  as 


(4.20) 

We  can  use  Equation  (4.20)  to  write  the  phase  velocity  for 
this  single  wave  as 


“k  +  A<uk 

R 


(4.21) 


in  terms  of  the  modulus  of  the  mode  amplitude.  From  Equation 
(4.17)  it  is  clear  that  ck  =  qk+)  (t«0)  ,  so  that, 

lck  1 2  *  !<Ik+;  <t=0>  I2  =  k2  c2  (4. 


where  we  have  used  Equation  (3.16)  to  write  the  second  equality 
in  equation  (4.22)  .  The  phase  velocity  may  therefore  be  writ¬ 
ten  as 

cph  "  j1  +  7  k2  <2j  (4.23) 

which  is  seen  to  agree  to  second  order  with  the  phase  velocity 
of  a  Stokes  wave,  i.e.,  have  the  appropriate  amplitude  dependence. 
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5.  SEA  SURFACE  CONTOUR  AND  ENERGY 

Consistent  vith  the  spirit  of  our  technique,  which  assumes 
the  non-linear  couplings  to  be  relatively  weak,  we  shall  cal¬ 


culate  the  wave  energy  to  only  second  order  in  the  wave  ampli¬ 
tude.  The  kinetic  and  potential  energies  per  unit  surface 
area  are,  respectively, 


and 


and 


K.E. 


pu  dz  , 


v  5 . 1 ) 


P.E. 


c/dx/ 1 

o  T  -<® 

1  f  1  2 

Lj  j  pg  C  dx  .  (5.2) 

o 

Use  of  Equations  (3.10)  and  (3.11)  gives  us,  to  second  order. 


K.E.  ss 


P.E. 


(5.4: 


Here  we  have  also  made  the  approximation  of  replacing  Gk  by  g  in 
Equation  (3.10).  The  distinction  between  Gk  and  g  is  important 
only  for  the  shorter  wavelengths,  which  are  expected  to  con¬ 
tribute  relatively  little  to  the  energy  and  also  to  violate  the 
expansion  in  Equation  (3.17). 


for  right- travel ling  waves  only,  the  energy  is  just 


K.E.  +  P.E. 


(5.5) 


V 


I 

4 


i 


< 
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The  surface  contour  obtained  from  the  definition  of  the 


mode  expansion  [Equation  (2.14)],  is, 


Ux,t) 


Im 


i [kx  -  w^t] 


(5.6) 


for  waves  travelling  only  to  the  right,  we  keep  only  terms  with 
k >  0,  here.  For  the  case  of  a  single  mode,  Equation  (5.6)  re¬ 
duces  to  Equation  (3.16),  of  course. 


We  wish  to  write  Equation  (5.6)  so  as  to  remove  the  pri¬ 
mary  wave  component  from  the  oscillation.  To  do  this  we  write 

the  primary  wave  number  k-k  so  that, 

o 

l*>k  *  +  IC )  M  (i)j^  |1  +  y  jj—  +  .  .  . 


where  the  group  velocity  (vG)  of  the  primary  wave  is  given  by 
7  The  wave  envelope  can,  therefore,  be  written  as. 


C (x,t) 


(+}  iKlx_vGtl 

e  e 


(5.8) 


Alternatively  we  can  express  Equation  (5.8)  in  terms  of  a  phase 
and  amplitude  modulation  by  defining  the  quantities 
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CR  +  1  GI 


,  .  itc  (x  -  v_t) 
e 


(5.9) 


and 


k  x  -  a)  t 
o  o 


(5.10) 


Substituting  Equations  (5.9)  and  (5.10)  into  (5.8)  yields 


C(x,t)  "  +  G?  j-^rcose  + 

▼  R  I  I  I G |  o  |Gl 


or 


| G |  Sin  6o 


C(x,t)  -  | G I  cos 


0^  -  arctan 
o 


ft)]- 


(5.11) 


(5.12) 


The  modulation  of  the  surface  waves,  i.e.  ,  the  surface  envelope 
is  therefore  given  by  the  function, 


I G  i 


/+\  i<(x-v  t) 

«k  *  e 


(5.13) 
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6.  SAMPLE  CALCULATION  AND  CONCLUSION 

The  purpose  of  doing  a  sample  calculation  is  twofold: 

(x)  to  test  the  numerical  methods  in  a  simple  case  which  can 
be  compared  to  an  analytic  calculation,  and  (ii)  to  deter¬ 
mine  instructive  ways  of  presenting  the  results  of  the  cal¬ 
culation.  For  our  test  calculation  we  have  selected  what 
we  feel  to  be  the  simplest  problem  which  is  still  of  some 
physical  interest.  The  initial  state  in  the  test  problem 
consists  of  three  equally  spaced  non-zero  modes,  the  center 
mode  being  an  order  of  magnitude  larger  in  amplitude  than 
its  side  bands.  The  interaction  between  the  modes  is  de¬ 
scribed  by  Equation  (3.24)  and  the  growth  of  all  the  modes 
in  the  system  is  calculated. 

This  problem  is  the  one  considered  both  experimentally 
and  theoretically  (linear)  by  Benjamin  and  Feir  (1967)  and 
which  lead  them  to  the  conclusion  that  a  Stokes  wave  is  un¬ 
stable.  The  experiment  consisted  of  generating  a  mechanical 
wave  of  fairly  large  slope  (0.17)  in  a  wave  tank  and  modula¬ 
ting  this  wave  with  two  low  amplitude  waves  (perturbations) 
at  the  side  bands  of  the  primary  wavenumber.  These  perturba¬ 
tions  were  found  to  grow  exponentially  from  out  of  the  back¬ 
ground  noise  on  the  tank  surface  and  eventually  caused  the 
primary  wave  to  break  up.  The  details  of  this  break-up  pro¬ 
cess  will  be  discussed  below  in  terms  of  the  present  calcula¬ 
tion. 
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We  have  limited  the  total  number  of  modes  in  this  first 
calculation  to  nine.  This  number  is  somewhat  arbitrary  but 
it  is  large  enough  to  indicate  the  general  behavior  of  the 
system  when  the  number  of  modes  is  increased.  Comparisons 
with  the  results  from  a  second  calculation  with  thirteen  modes 
is  also  made.  Also,  the  present  number  of  modes  does  yield 
accurate  quantitative  results  over  a  time  scale  long  compared 
to  the  characteristic  growth  time  of  each  of  the  separate  modes. 

The  wavenumbers  for  our  modes  are  chosen  such  that, 

k  *  17  <  i  K  -  0,  1,  2,  .  .  .  ,  N  (6.1) 

where  L  is  a  length  representing  the  region  of  interest  of  the 

ocean  surface  and  k  is  an  integer.  Because  we  wish  to  make 

some  comments  on  the  comparison  of  the  present  calculation  and 

the  experiment  of  Benjamin  and  Feir,  we  select  a  length  (L  = 

4 

2.1946  x  10  cm)  and  integer  (<  ■  100)  such  that  the  central  wave- 
number  (k)  is  that  of  the  primary  wave  in  their  experiment, 
i.e.,  k*  0.02863  cm  Our  first  calculation  will  not  use  the 
experimental  amplitudes  since  these  values  lead  to  a  very 
rapid  break  up  of  the  primary  wave.  The  initial  configuration 
is  shown  in  Figure  3  labeled  "time  «  0  sec."  We  have  selected 
an  initial  amplitude  of  the  primary  mode  to  be  |q^+)  |  =0.071  and 
I  for  the  secondary  modes  which  corresponds  to  wave 

amplitudes  of  2.46  cm  for  the  primary  and  0.45  cm  for  the  greater 
and  0.55  cm  for  the  lesser  of  the  secondary  modes.  The  remainder 
of  the  nine  node  amplitudes  are  initially  zero. 

I 
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The  results  of  this  calculation  are  presented  in  Figures 
1  through  6.  In  Figure  3,  we  show  a  series  of  snapshots  of 
the  modulus  of  the  mode  slopes.  Each  snapshot  indicates  how 
the  energy  has  redistributed  itself  between  the  modes  from  the 
preceding  snapshot.  We  can  see  that  the  energy  of  the  system 
which  at  time  zero  was  concentrated  in  the  three  central  modes 
diffuses  outward  into  the  neighboring  modes  of  the  system  in 
time.  Those  modes  closer  to  the  central  mode  grow  faster  than 
the  more  distant  modes.  It  is  evident  from  Figure  3  that  the 
evolution  of  the  modes  closely  resembles  the  process  of  dif¬ 
fusion.  This  notion  of  the  diffusion  of  energy  between  modes 
due  to  non-linear  interactions  will  be  explored  more  fully  in 
a  subsequent  report  (PD-72-029)  which  models  the  interaction 
process  between  ocean  waves  in  terms  of  a  transport  equation. 

In  Figure  4  the  sensitivity  of  the  preceding  calculation 
to  changes  in  mode  number  is  shown.  Maintaining  the  spacing 
between  modes  we  increase  the  width  of  our  k-space  interval 
by  four  modes  maintaining  the  same  central  mode.  It  is  clear 
that  these  additional  modes  have  no  effect  on  the  growth  for 

t  <  80  sec  since  only  the  12th  mode  has  significant  amplitude 
up  to  this  time. 

Figure  5  is  a  continuous  representation  of  the  informa¬ 
tion  in  Figure  3,  showing  how  each  of  the  mode  amplitudes 
grow  as  a  function  of  time.  We  can  see  from  Figure  3  that 
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FIGURE  (3) 


FIGURE  m 


Wtf  =  lM0DE  SLOPEl2 


the  side  bands  do  not  grow  at  exactly  the  same  rate,  but  close 
enough  that  a  single  line  may  be  used  in  Figure  5  to  represent 
the  growth  of  both  these  modes.  The  central  mode  (kq  =  100)  is 
labeled  "Primary  Mode"  in  Figure  5,  the  first  side  bands  (k  =  99 
and  101)  are  labeled  "Secondary  Mode,"  etc.  In  this  figure  we 
compare  our  calculations  to  an  analytic  result,  i.e.,  the  Born 
approximation.  The  Born  approximation  for  this  problem  is 
calculated  by  assuming,  (i)  that  the  amplitude  of  the  primary 
mode  is  constant  throughout  the  interaction  process,  and  (ii) 
the  frequencies  of  the  primary  and  secondary  modes  are  approx¬ 
imately  equal.  The  interaction  equations  then  become, 

i^'x  p 

%+l  ‘JKo+1+%  <0-l]  '  (6.2) 


and 


q  a  0 

^  ir 

Ko 


q  .  a 
k  -1 


1(jV 


■Vl+q 


K0 


since,  by  assumption  (ii) 


(6.3) 

(6.4) 


a  w 


Ko+1 


(6.5) 


We  may  use  Equations  (6.2)  -  (6.4)  to  construct  the  equation 

3t  'Vi!2'  ft  'Vi!''  -  “toqt0I",{qK0+l‘3K0-l)  (6-6) 

where  Im|  }  refers  to  the  imaginary  part  of  the  bracketed  quan¬ 
tity.  Since  our  initial  conditions  are  q  (0)  =  q  (0)  =  q 

Ko~1 

and  q^  (0)  =  q*  (0),  Equation  (6.6)  has  the  solution 
o  o 
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%+i(t)  ■  qK0-i!t)  ■  qHhh0K/l]  ■  isinhh<0K0|2£]j 

(6.7; 

The  Born  approximation  to  the  secondary  modes  plotted  in  Fig¬ 
ure  5  is,  therefore, 

+J2  =  -1^  =  q2  cosh[  “k  I<3k  I2fcl  •  (6*8) 

We  see  that  the  growth  of  the  secondary  mode  is  similar 
to  Equation  (6.8)  in  the  region  where  the  approximations  made 
are  nearly  valid,  that  is,  for  a  near  constant  primary  mode. 

The  development  is  markedly  different,  however,  after  the  higher 
modes  have  grown  to  an  appreciable  size,  that  is,  the  secondary 
modes  cease  growing  and  begin  losing  their  energy.  The  energy 
is  depleted  from  the  secondary  and  primary  modes,  causing  the 
growth  of  the  more  central  modes  to  be  inhibited.  It  is  not 
shown  here,  but  the  energy  drain  of  the  primary  mode  does  not 
persist,  it  begins  growing  after  the  other  modes  of  the  sys¬ 
tem  are  of  the  same  order  of  magnitude.  The  detailed  manner 
of  this  growth  hes  yet  to  be  explored. 

In  Figure  6  we  illustrate  the  envelope  of  the  surface 
waves  in  a  region  of  ocean  438  meters  long.  The  envelope  is 
calculated  as  it  travels  with  the  group  velocity  of  the  pri¬ 
mary  wave  (9.26  m/sec)  so  that  each  envelope  shown,  i.e.,  dif¬ 
ferent  time  snapshots,  is  constructed  from  the  same  group  of 
waves  as  they  propa3ate  along  the  ocean  surface.  In  Figure  6a 
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AMPLITUDE  OF  SURFACE  ENVELOPE 


(6a) 


(6b) 


DISTANCE  IN  UNITS  OF  Xprjmary  (=  2.194  meters) 

(6c) 

FIGURE  (6) 
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we  see  that  (time  t=0)  the  modulation  of  the  ocean  surface 
is  slow,  but  significant  over  the  region  of  interest.  This  is 
due  to  the  large  initial  amplitude  of  the  secondary  modes. 

The  step  size  is  given  in  units  of  the  primary  wavelength,  i.e., 
1  =  2.194  meters.  The  initial  undulation  of  the  envelope  seems 
to  be  compressed  in  time  into  a  series  of  bumps.  These  bumps 
rapidly  become  accentuated,  forming  wave  packets  on  the  sur¬ 
face.  It  will  be  seen,  however,  that  the  large  structure  in 
Figure  6  is  due  to  the  periodic  boundary  conditions  imposed 
on  the  problem  and  does  not  represent  the  breakup  seen  ex¬ 
perimentally  by  Benjamin  and  Feir. 

The  significant  structure  in  Figure  6  lies  between  the 
modulation  peaks.  In  this  region  the  modulation  remains 
fairly  constant  during  the  time  in  which  energy  is  diffusing 
out  of  the  center  mode  (see  Figures  3  and  4).  As  the  c.mpli- 
tude  of  the  higher  modes  increases,  however,  the  structure  of 
the  wave  envelope  changes.  The  extreme  case  is  shown  in  Fig¬ 
ure  6c  where  the  detailed  modulation  of  the  surface  waves 
has  clearly  developed.  It  is  difficult  to  determine  the 
surface  structure  by  looking  at  only  a  picture  of  the  modula¬ 
tion  so  we  introduce  a  parameter  to  indicate  the  degree  of 
distortion. 

To  characterize  the  distortion  of  the  ocean  surface,  we 
use  the  ratio  of  the  difference  between  the  maximum  and  mini- 
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mum  envelope  height  to  the  average  envelope  height.  In  Fig¬ 
ure  7  we  use  this  quantity  to  indicate  the  growth  in  the  sur¬ 
face  wave  amplitude  modulation.  The  solid  curve  refers  to 
the  sample  calculation  above,  the  dashed  to  a  simulation  of 
the  Benjamin-Feir  experiment.  We  see  in  Figure  7  that  for 
t^70  sec,  essentially  no  change  is  observed  in  the  modula¬ 
tion  of  the  surface.  Note  that  we  have  removed  the  effects 
of  the  large  bumps  in  Figure  6  and  are  concerned  only  with 
the  "actual"  modulation  which  lies  between  the  bumps.  The 
sudden  growth  of  the  surface  distortion  (t>70  sec)  occurs 
when  the  magnitude  of  the  modes  becomes  comparable  (see  Fig¬ 
ures  3  and  4)  so  that  large  transfers  of  energy  can  take 
place  in  small  interaction  times. 

The  second  curve  in  Figure  7  represents  the  growth  of 
the  surface  distortion  from  the  calculation  using  the  ini¬ 
tial  conditions  of  the  B-F  experiment.  This  calculation 
again  uses  nine  modes,  with  the  central  mode  corresponding 
to  the  mechanically  generated  wave.  The  experimental  condi¬ 
tions  were  simulated  by  giving  each  of  the  remaining  modes 
the  uniformly  small  value  o.  0.001,  which  is  intended  to  re¬ 
produce  the  noise  at  the  tank  surface.  In  terms  of  our  vari¬ 
able,  the  central  mode  amplitude  is  q^+p  =  0.17  corresponding 
to  a  primary  wave  of  5.93  cm  amplitude.  These  initial  con¬ 
ditions  are  shown  by  the  "time  =  0"  graph  in  Figure  8. 
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FIGURE  (7) 
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TIME  IN  SECONDS 


FIGURE  (8) 
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Figure  8  depicts  the  modulus  of  the  mode  amplitude  at 
different  instants  of  time  just  as  Figure  3  for  the  test  cal¬ 
culation.  e  can  see  that  the  non-linear  interactions  pref¬ 
erentially  amplify  the  side  bands  which  differ  from  tne  pri¬ 
mary  by  ±3Ak,  where  Ak  is  the  step  size  in  wavenumber  space. 
These  side  bands  are  1.3  and  0.7  of  the  central  wavenumber 
(iCq)  •  This  is  in  essential  agreement  with  the  perturbation 
analysis  of  Benjamin  which  shows  that  the  frequency  side  bands 
at  1.1  and  0.9  of  would  be  preferentially  amplified  from 

out  of  the  background  noise.  The  effect  of  increasing  the 

number  of  modes  in  this  calculation  is  shown  in  Figure  9, 

where  we  have  reduced  the  step  sizr  Ak  by  a  factor  of  two. 

In  Figure  9  the  mode  numbers  which  differ  from  the  primary 

by  ±10%  are  modes  1  and  13.  We  can  see  that  instead  of  a 

single  mode  on  each  side  of  the  primary  being  preferentially 

amplified  we  have  a  preferred  group  of  modes  being  amplified. 

This  would  seem  to  indicate  the  validity  of  the  preferred 

mode  concept  in  the  continous  limit,  where  a  group  of  waves 

in  a  region  6<u  about  the  average  position  ±10%  of  to  would 

o 

be  picked  out  of  the  background. 

In  Figure  10  we  again  illustrate  the  modulat-on  of  the 
surface  envelope  on  a  stretch  of  ocean  200  times  the  primary 
wavelength  (Xp*  2.194  meters).  Comparing  the  initial  distor¬ 
tion  with  that  of  the  test  calculation  in  Figure  6a,  we  see 
that  the  direct  comparison  is  somewhat  deceptive  in  that  the 
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B“F  case  has  a  greater  initial  distortion  which  is  not  evident 
from  the  figure.  In  Figure  10  the  surface  modulation  for  both 
the  9  and  13  mode  B-F  cases  discussed  are  given.  We  see  that 
the  effect  of  decreasing  the  mesh  size  by  a  factor  of  two  in 
k-space  has  displaced  the  central  bump  in  the  9  mode  calcula¬ 
tion  from  100  Ap  to  20C  Ap  in  the  13  mode  calculation.  The 
effect  of  this  reduction  is  most  clearly  seen  in  the  "time  = 

30  graph.  Although  the  9  mode  calculation  does  show  the 
modulation  of  the  surface  waves,  the  envelope  becomes  distorted 
as  it  approaches  the  1/Ak  =  100Ap  point  in  the  figure,  i.e., 
the  periodic  boundary.  This  effect  is  markedly  reduced  in  the 
13  mode  calculation  which  shows  precisely  the  type  of  surface 
modulation  observed  by  Benjamin  and  Feir.  This  comparison 
indicates  that  caution  must  be  exercised  in  the  selection  of 
the  mode  spacing  used  in  a  particular  calculation.  The  am¬ 
plitude  scales  in  Figure  10  are  arbitrary  because  the  9  and 
13  mode  curves  were  shifted  so  as  to  provide  the  best  visual 
comparison.  The  initial  curves  are  actually  superimposed. 

The  distance  the  primary  wave  travelled  in  the  B-F  ex¬ 
periment  prior  to  breaking  up  was  approximately  twenty-eight 
wavelengths.  The  corresponding  time  interval  is  33  seconds. 

We  see  from  Figure  10  that  the  distortion  of  surface  becomes 
quite  significant  (a  10%)  at  this  time,  which  agrees  well  with 
the  experiment. 
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It  is  clear  that  when  the  primary  wave  initially  contains 
most  of  the  energy  in  the  wave  system,  that  the  coherent  non¬ 
linear  interactions  between  the  surface  gravity  waves  lead  to 
a  breakup  of  this  wave  into  packets.  The  characteristic  time 
for  this  breakup  is  determined  by  the  details  of  the  initial 
conditions.  The  breakup  is  seen  as  a  distortion  of  the  sur¬ 
face  envelope.  This  distorting  occurs  as  energy  coherently 
diffuses  from  the  primary  to  the  other  modes  in  the  system. 

As  additional  modes  develop  the  surface  distortion  increases 
markedly,  but  the  pattern  of  the  distortion  does  not  change. 

We  see  this  in  Figure  10  where  the  initial  pattern  in  the  B-F 
experiment  translates  and  grows  but  does  not  change  in  shape. 
This  would  indicate  that  the  group  velocity  of  the  pattern, 
which  is  that  of  the  primary  wave,  changes  only  slightly  as 
the  other  modes  develop;  also  that  the  phases  of  the  other 
modes  "lock"  onto  that  of  the  primary  so  as  to  create  the 
stationary  pattern. 

If  such  an  effect  were  present  on  the  real  ocean  it 
could  be  observed  with  radar.  The  radar  return  from  such  a 
surface  would  indicate  the  envelope  structure  and  "see"  the 
phase-locking  effect.  The  distorted  envelope  would  yield  a 
strong  correlation  centered  at  the  frequency  of  the  primary 
wave. 


60 


APPENDIX  A 


To  integrate  Eq.  (3.29)  over  the  interval  0<t<T, 
we  choose  a  time  step  A,  where  T/A  is  an  integer.  Let 


tj  =  jA ,  j=0,  1,  2,  ...  ,  T/A  . 


Integration  of  Eq.  (3.29)  over  the  interval  t^  to  t^+^  is 
done  in  the  form: 


q^+)  (t.  +  A)  = 


i  ...  /  aiM  .  VA 

\  2  “k  / 


^+)(V 


T, 

+  i  * .  .  _  r 


£+p=k+n  £pn;k  £pn;k 


(A.l) 


where 


I)lpn;k(:i) 


f 

t . 
3 


fcj+A 


H.  e 
S,pn 


^Jlpnk11 


dt 


(A. 2) 


The  numerical  evaluation  of  Eq.  (A.  2)  is  done  as 
follows : 

Write  [we  drop  for  the  moment  the  (&,p,n,k)  subscripts] 


H  ( t)  s  Hm(t)  e1<Mt) 


(A. 3) 


where  and  $  are  real.  We  suppose  for  the  evaluation  of  I(j) 
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the  q,[+^  (t . )  and  cr/+^ 


'fcj'  an<*  ^k  'fcj  ”  A  )  are  stored.  We  can  thus  evaluate 


wf  ) 


sW+fHi  ' 


tH-a(tj>  -  Hm(tj  -  A)]/A  ' 


(A. 4) 


-Kt.  +  j)  s  4>(t.)  +  | 


>!  -  l$(t.)  -  *{tj-A)]/A 


(A. 5) 


We  may  then  express  I(j)  [Eq.  (A. 2)] 


I  ( j  )  =  e 


it« (t j  +  |  )+*Q] 


i[Q+4>  ]x 

[Hq  +  H1t]  e  1  dx 


i[fl(t.  +  *)  +  *  ] 

1  2  O 


=  e 


'Vo  +  H1L1> 


where 


=  TTOTT  sin  C(n  +  VA/2]  » 

=  - - - 5-  \  2  sin  [(a  +  4> , )  A/2] 

(n+^r  (  1 

-  ((fi  +  G^a)  cos  l(n  +  *x)  A/2]  | 


(A. 6} 


(A. 7) 


In  Equation  (A. 6)  we  have  introduced  the  new  variable 
of  integration  x  =  t  -  (tj  +  A/2). 

To  check  the  numerical  accuracy  of  the  calculation  at 
each  step  of  the  integration  the  energy  is  calculated  using 
Equation  (5.5).  Since  the  energy  will  be  conserved  in  the 
absence  of  wind  and  viscosity,  this  provides  a  quantitative 
check  of  our  procedure. 


61 


REFERENCES 


Benjamin,  T.  B. ,  (1967),  Royal  Soc.  of  London,  Proceedings , 

A299,  59. 

Benjamin,  T.B-  and  J.  E.  Feir,  (1967),  J.  Fluid  Mech.  2J_,  417. 

Cohen,  B.  I.,  A.  N.  Kaufman  and  K.  M.  Watson,  (1971),  Phys . 

Rev.  Lett.  29,  581. 

Hasselmann,  K. ,  (1961),  J.  Fluid  Mech.  12,  481. 

Hasselmann ,  K. ,  (1966),  Reviews  of  Geophysics  4 ,  1 . 

Phillips,  0.  M. ,  (1966),  The  Dynamics  c.  the  Upper  Ocean , 

Cambridge  University  Press. 

Whitham,  G.  B.,  (1966),  Royal  Soc.  of  London,  Proceedings ,  V299 . 


65 


